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ABSTRACT: We describe a technique for constructing a symplectic transfer map for a charged par- 
ticle moving through an accelerator component with arbitrary three-dimensional static magnetic 
field. The transfer map is constructed by symplectic integration; by representing the map at each 
step of the integration by a mixed-variable generating function, exact symplecticity is ensured. By 
using an appropriate integration algorithm, there is no necessity to make the paraxial approxima- 
tion. The technique is illustrated by application (in one degree of freedom) to a quadrupole magnet 
with strong octupole component and fringe field. 
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1. Introduction 

The ability to compute accurately and efficiently the trajectory of a particle in an accelerator beam 
line is fundamental to accelerator design and operation. Where the beam line includes components 
with complex electromagnetic fields, and where fast tracking and a high degree of accuracy are 
required, this problem presents significant challenges. Considerable work over the years has led to 
the development of a variety of techniques for particle tracking (see, for example, |jl|-||]); however, 
since the equations describing particle trajectories cannot be solved exactly, various approximations 
are always needed. Often, there is a compromise that must be made between accuracy and speed. 

There are two approaches to tracking a particle through the electromagnetic field in a given 
accelerator component: either, the equations of motion may be integrated numerically, or, a transfer 
map may be applied for the entire component. Although the equations of motion are well known, 
numerical integration is in general too slow for long-term particle tracking in a storage ring. A 
transfer map, consisting of a set of equations relating the particle position and momentum at the 
exit of the component to the position and momentum at the entrance, provides a much more efficient 
tracking method; however, the problem then is to construct with good accuracy the transfer map for 
a given component. One approach is to integrate the equations of motion using a differential algebra 
(DA) code to produce the transfer map in the form of a multi- variable power series. Although 
a map in this form is very convenient for particle tracking, in general the power series must be 
truncated at some order in each step of the integration, to keep the number of terms manageable. 
As a consequence of the truncation, the final map will not be symplectic: this may have an impact 
on detailed aspects of the particle dynamics (for example, whether or not a particle remains in 
a storage ring over a large number of turns Techniques have been devised to construct 
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symplectic, finite-order Taylor maps (Cremona maps) that correspond to a symplectic Taylor map 



that has been truncated at some order [hQ, 11]. However, these methods are not especially easy to 



implement, and the consequences of modifying a map to make it symplectic are not always clear 

The need for efficient symplectic tracking tools has motivated the development of techniques 
for manipulating symplectic transfer maps (symplectomorphisms) using, for example, representa- 
tions commonly known as Lie transformations [Q, [12|]. Venturini and Dragt have shown how to 
construct the transfer map for an accelerator component in the form of a Lie transformation, using 



detailed magnetic field data Q13|]. While Lie transformations provide a valuable analytical frame- 
work, it is not straightforward to implement the relevant formulae in an accelerator simulation 
code. As an alternative, we discuss here the construction of a mixed-variable generating func- 
tion in particular, to represent the dynamics in an accelerator component with some specified 
static inhomogeneous magnetic field. Mixed-variable generating functions are already used for 
particle tracking where symplecticity is important [|7j |9|, [l^- 21 ]. Besides ensuring symplecticity, 
mixed-variable generating functions also have potential advantages over Taylor maps, in reducing 
the number of coefficients needed to describe the map, and in improving the accuracy of a map 
for a given accelerator configuration that is obtained by interpolation between known maps rep- 



resenting a set of "reference" configurations [ |22| , |23[] . However, the common approach to using 
mixed- variable generating functions is first to construct truncated multi- variable Taylor expansions 
for the canonical variables (approximating the transfer map, usually for an entire turn of a storage 
ring) and then to construct the generating function that gives this map to some desired order. Such 
a process can be used to "symplectify" a truncated power series; but then the same concerns ap- 
ply as in the case of construction of Cremona maps from truncated power series. A more direct 
approach, that we consider here, consists of working entirely ab-initio with mixed-variable gener- 
ating functions. The advantage is an improvement in efficiency and accuracy that can be achieved 
by avoiding conversions between different (and, in some cases, non-symplectic) representations of 
the same map. 

Although the manipulation of maps in the form of mixed-variable generating functions is not 
trivial, it is at least possible to obtain from the generating function an exact set of algebraic equa- 
tions relating the initial and final values of the dynamical variables that preserves the symplectic 
structure. Unfortunately, these equations then need to be solved numerically for given initial con- 
ditions, but this can be done to any desired precision with reasonable efficiency, using an iterative 
technique such as Newton's method. The symplecticity of the map can then be preserved, to ma- 
chine precision. 

The above arguments suggest that a technique for constructing the generating function for a 
given accelerator component with specified static magnetic field offers a potentially useful tool in 
particle tracking: we describe such a technique in Section ||. We see that relatively straightforward 
expressions can be written down for the mixed-variable generating function for each integration 



step in the case that the paraxial approximation can be made for the Hamiltonian (Section 2.1); 



slightly more complicated (though still very manageable) expressions can be used where a more 



accurate form of the Hamiltonian is needed (Section 2.2). A key step in the integration is the 



composition of two generating functions: this is discussed in some detail in Sections ^31 and 2.4, 



where we present two methods for performing the composition. Then, in Section ^ we illustrate 
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the techniques developed in Section || by applying them to the case of a magnetic quadrupole with 
strong octupole component and fringe field. 



2. Symplectic integrators 

We shall consider the case of a charged particle moving through a static magnetic field B = V xA. 
Let s parameterise the reference trajectory, which by assumption is an orbit in the magnetic field. 
The canonical variables on the (extended) phase space are (x,s,p,p s ) where x = (x,y,z) and p = 
{pxiPyiPz)- x an d y denote the position of the particle in a plane perpendicular to the reference 
trajectory at s, which the particle crosses at time t; and z = s/po — ct where po is the speed (as 
a fraction of the speed of light) of a particle with a fixed reference (kinetic) momentum Pq. The 
transverse components of the canonical momentum are defined by: 

Y(\p\)mx + qA x (x,s) 

Px = — = , (2.1a) 

"o 

y(\P\)my + qAy(x,s) 
Py = — — " p , (2- lb) 

where y(|j3|) is the relativistic factor; the particle velocity is j8c = (x,y,s); m and q are the mass 
and charge of the particle; and the dot indicates a derivative with respect to time. The momentum 
conjugate to the longitudinal co-ordinate z is the energy deviation p z , defined by: 

E 1 no\ 

Pz = s r-> (2 - 2 ) 

Poc Po 

where E = y(|/$ |)mc 2 is the kinetic energy of the particle. The reference momentum can be written 
as Po = fioYomc, where yo = y(|po|). Note that third component of the particle co-ordinate vector 
x = (x,y,z) actually specifies the time at which a particle arrives at a given longitudinal position 
within the magnetic field; so for a magnetostatic field, A is independent of z. 
It is convenient to work with the normalised vector potential, defined by: 

a = A (2.3) 
Bp 

where Bp =P /q is the beam rigidity and the vector a has Frenet components (a x ,a y ,a s ). 

In these variables, the Hamiltonian for particle motion in a static magnetic field with a straight 
reference trajectory is: 



(2.4) 



H=P i~iik +p ) ~ {Px ~ ax)2 ~ {Py ~ ay)2 ~w ~ as + Pi ' 

where y> = y(po). The independent variable is denoted by a. Note that: 

ds dH _ 1 
da dp s 

and therefore: 

s = o + o , (2.6) 
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where Ob is an (arbitrary) constant of integration. 

The method will be illustrated for a particle that moves with y = p y = p z = 0, i.e. we consider 
only the evolution of the variables x,p x and s. The generalisation to the full number of degrees of 
freedom is straightforward, and requires no significant new features. In addition, we assume that 
the particle is ultra-relativistic, so that we can take the limit /o — > °°. By making an appropriate 
gauge transformation, it is always possible to work in a gauge in which one component of the vector 
potential vanishes at all positions in space: we shall choose a y = 0. The Hamiltonian is then: 



H 



1 " (Px 



Our assumption that y = p y = as the particle moves through the field requires that: 



da x 
dy 



y=0 



da s 
~dy~ 



0. 



(2.7) 



(2.8) 



y=0 



This is a restriction on the field, but again the generalisation to include other cases requires no 
significant new features. 

An accelerator component is described by specifying the components of the vector potential. 
Neglecting radiation and collective effects, the transfer map for a particle in a magnetostatic field is 
symplectic. Therefore, the values of the canonical variables at the exit of the accelerator component 
can be related, at least in principle, to the values of the variables at the entrance by a mixed- variable 



generating function. Following Goldstein's nomenclature Q14[], with (x\,p x \) the entrance variables 
and {x2,p x i) the exit variables, we write F{x\,p&) as a generating function of the second kind. 
Then, the entrance and exit variables are related by: 

dF 



X2 



Pxl 



dpxl 

dF 

dx\ ' 



(2.9a) 
(2.9b) 



Our objective is to develop a procedure for constructing F{x\,p x 2) in a given static magnetic field. 
2.1 Integrator for the Hamiltonian in the paraxial approximation 

Assuming that the dynamical variables and normalised vector potential have small values, we can 



make the paraxial approximation for the Hamiltonian given by Eq. (2.7) 

f-a s + p s . 



1 



Hk, -\ + -( Pr -a r Y-a, + p*. (2.10) 
In this case, a method for constructing a convenient symplectic integrator has been described by 



Wu, Forest and Robin [24]. A more general symplectic integrator, not requiring the paraxial ap- 



proximation, is described in Section 2.2. However, where it is possible to make the paraxial approx- 
imation without losing the desired level of accuracy, the advantage of doing so is that integration 



can be performed more simply. Following [24], using the Hamiltonian in the paraxial approxima- 



tion Eq. (2.10), the integration of the particle motion through a step A a (corresponding to a step 
As = Act in the field) can be expressed as a sequence of transformations, for example: 



^# = o ^2 -^3 -^2 1 ° ^\ 



(2.11) 



-4- 



where: 



^is = s + ^Ao, (2.12a) 

Jift = ^i^AcT^, (2.12b) 
2 

^2Av = Px~a x , (2.12c) 

•^f 1 Px = Px + ax, (2.12d) 

^# 3 x = x + Aop x . (2.12e) 

The effect of the transformations ^#2, -^3, on variables not explicitly shown is to leave them 
unchanged. These transformations provide a second-order integrator (the error is 

If the components of the vector potential a are expressed as power series in x and s, then at 
this stage we can use a differential algebra code to compose the transformations first for each step 
through the field, and then for successive steps through the field. However, except for trivial cases, 
this will require truncation of the power series in order maintain a manageable number of terms in 
the power series, and the truncation will mean a loss of symplecticity. Furthermore, we know from 
the symplectic condition that many of the coefficients in the power series will be related to each 
other: this direct approach is therefore more computationally expensive than necessary. 

Instead, we can write down mixed-variable generating functions corresponding to each of the 
individual transformations. Using generating functions of the second kind, it can be seen that 
generators for the above transformations may be written: 

Fi{x\,Px2iSi,p s2 ) = xip X 2-^Aca s (xi,si) + sip s2 + }^Acp S 2, (2.13a) 

F2{x\,p x2 ,si,p s2 ) = xip x2 + / a x (x,si)dx + sip s2 , (2.13b) 

Jo 

^ixuPx2,s\,p S 2) = xip x2 + -Aap 2 x2 + sip s2 , (2.13c) 

where F, is the generator for the map Mi. The generator for M^ is: 

- 1 f Xl 
F 2 (xi,p x 2,si,p s2 ) =x\p X 2- I a x (x,s\)dx + s\p s 2. (2.14) 

Jo 

To proceed, we need a way to compose the generating functions. We discuss two possible methods 



for this in Sections g3| and [2.4l However, we first discuss an integrator, expressed in terms of mixed 
variable generating functions, that can be used even in cases where the paraxial approximation is 
not valid. 



2.2 Symplectic Runge-Kutta integrator 

In the previous section, we wrote down expressions for mixed-variable generating functions for 
integration steps where the dynamics could be described by a Hamiltonian in the paraxial approxi- 
mation. While the expressions obtained can be expressed rather simply, there are cases of interest 
where the paraxial approximation is not valid. For such cases, at the cost of rather more compli- 
cated formulae, we can use an alternative integration method, such as a symplectic Runge-Kutta 
algorithm. 
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For a Hamiltonian system, a symplectic Runge-Kutta method may be written [g5|] : 

dH{x\p x ) 



x 2 = xi+Aaj^bi- 

i=i 



Pxl = Pxl 



Aa£, 

7=1 



dp x 



(2.15a) 
(2.15b) 



where the "intermediate" variables (x',p' x ) are given by: 

dH{xi,pi) 



x' = x\ + Aa a,- 



Aa ^ cij 

7=1 



dH(xJ,p J x ) 



The transformations ( [2.1 5| ) and ( |2.16| ) are symplectic if: 

biaij+bjQji b,b h 



(2.16a) 
(2.16b) 

(2.17) 



for all i and /. A simple case is given by n = 1, an = 1/2 and b\ = \: this provides a second-order 
symplectic integrator. A fourth-order integrator is given [p6|] by « = 2, and: 



0*7 ) 
(*/) 




(2.18a) 
(2.18b) 



The transformations in Eq. (2.15) can be derived from a mixed- variable generating function of the 
second kind, given by [25]: 



A , , 7 A dH(xj,pi) dH(xJ,p£) 
F = x\ pxi + Aa £ ft,// (x ( , /£) - Aa 2 £ 6,^ — 



(2.19) 



i=l i,y=l 

To write the generating function explicitly, it is necessary first to solve for the intermediate vari- 
ables (x',p' x ) using Eq. ( 2.16| ), then use Eq. ( 2.15b ) to eliminate p x \ (thus expressing the generating 
function F purely in terms of x\ and Pxi). 

Although each integration step requires the solution of a number of equations in several vari- 
ables, it turns out that the solution can be accomplished in a straightforward fashion (using an itera- 
tive technique) with an algebraic code (such as Mathematica), or a DA code. We have implemented 
the fourth-order Runge-Kutta integrator (with coefficients ( 2.18 )), using the mixed- variable gener- 
ating function ( 2.19 ) in Mathematica; we return to the question of solving the systems of equations 
necessary to construct the mixed- variable generating function (in particular, Eqs. (2.16)) in Section 



2.4 



Note that Wu, Forest, and Robin [24] discuss the extension of their method for the paraxial 
case, to cases where the paraxial approximation is not valid. However, this requires (local) time 
integration, which adds complexity to the problem. The advantage of the Runge-Kutta integrator is 
that one can adopt a unified approach, using (effectively) path length as the independent variable, 
for Hamiltonians of sufficiently general form for most particle tracking applications in accelerators. 
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2.3 Generating function composition: induction method 

Let Fa(x\,Pjq) and Fb(x2,P&) be generating functions of the second representing the transforma- 
tions from (xi,p x i) to (x2,Px2) and from (x2,Pxi) to {xj,,Pxi) respectively. Our goal is to find a 
single mixed variable generating function Fc(xi,p x s) that represents the effect of the combined 



transformation of Fa and Fg. We first note that Fa(x\,p&) an d Fb(x2,P x j) satisfy fll4|]: 

cLFa 

3— = Px\X\ +X 2 px2, (2.20) 

da 

and 

3— = /?*2*2 + *3/?.t3- (2-21) 

where the dot now signifies differentiation with respect to the independent variable, a. If we define: 

F C = F A - x 2 p x2 + F B , (2.22) 

we see that Fc satisfies: 

= Pxixi + Xip X 3 ■ (2.23) 

do 

Therefore, if we express Fc as a function of x\ and p X 3, then Fc will provide a mixed variable 

generating function of the second kind, for a transformation from (xi,p x i) to {x$,pxi). In general, 



Eq. (2.22) gives Fc as a function of x\, p X 2, X2, and p x $. However we can use the equations: 



X2 = ^~, (2.24) 
aPxl 

and 

Px2 = 3—, (2.25) 

to eliminate and ^2. We are then left with Fc as a function of x\ and p V 3 only, as required for a 
mixed-variable generating function of the second kind. If Fa and Fb are written as power series in 
their respective variables, then elimination of p X 2 and X2 may be performed as follows. Equations 
( 2.24 ) and ( 2.25| ) may be written in the form: 



I &1(X2,PX2)=0, 

\&2(x 2 ,Px2)=0, ^ ■ ' 

where &\ and ^2 are polynomials of order N in the given arguments. In general it is difficult to 
find exact solutions to ( 2.26| ). However, we only require perturbative solutions to some order in a 



perturbation parameter e, where the dynamical variables are of order s. Specifically, e = implies 
(xi,p x i) = (0)0) an d we are on the ideal orbit. Thus we propose solving the perturbation equation: 



> l{x2,Px2) = O{S N+l 
&2(X 2 ,PX2)=0(S N+1 ) 



N + u (2-27) 



We define the vector X: 




X = z . (2.28) 
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Then we can re- write Eq. (2.27) as: 




(2.29) 



where each term in the summation is a vector, each component of which is a monomial of order r 
in the components of its vector arguments. Each of the T r can be regarded as a tensor of rank r + 1 , 
so that To is a 2-component vector, Ti is a 2x2 matrix, and so on. Using an index notation, we 
would write, for example: 



[f 3 (x 1 ,x 2 ,x 3 )].= £ [f 3 ] [x,] j [x 2 ) k [x 3 ] r 



(2.30) 



The T r are obtained from 3?\ and <^ 2 , and are symmetric (i.e. [T r (Xi , . . . ,X r )] . is unchanged if 
any two of the arguments are interchanged). 

We make two assumptions: first, that lim e _>o(Ti) is invertible; and second, that To = 0(e). 
Using the first assumption, we can define: 



T k = -t^t k . 
In Appendix |A|, it is shown that a solution to Eq. ( 2.29 ) is: 

N 

X=£c r , 

r=l 

where ci = To, and for r > 2: 



(2.31) 



(2.32) 



c r — T^(c m | , . . . , c mk ) . (2.33) 

k—2 m[,m2,....nif;—\ 
mH \-mk=r 

Here the inner summation is performed over all combinations (m\,m%, . . . ,rrik) where each m ; > 1 
and m\-\ Vm^ = r. Thus ( [2.33 ) can be rewritten as a multiple sum: 



r r r—m\r—m\—mi r—m\ m/,-2 

c r = 52 ^ ^ 52 ■'■ 52 T/.(c mi , . . . ,C mk _ l ,C r _ m , mt] 

jfc=2ff!l=lff!2=l I»3=l mt_]=l 



) 



The c, are given inductively, so that: 



c 2 = T 2 (T ,T ), 

c 3 = T 3 (T ,To,To) + 2T 2 (To,T 2 (To,T )), 



(2.34) 



(2.35) 
(2.36) 



and so on. We note that c r = 0(e r ). 

The above formulae may be generalised to any number of degrees of freedom. Thus, this 
procedure leads, in principle, to explicit expressions for the coefficients in Fc (expressed as a power 
series) in terms of the coefficients in Fa and F#. 

However, applied to the composition of generating functions, the solution expressed in equa- 
tions (2.32) and (233) is not purely numerical, since the coefficients in the polynomial equations 
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( J2.26[ ) must be expressed in terms of dynamical variables (in the case of one degree of freedom, x\ 
and Pxi)- Unfortunately, for all but the simplest cases (low-order maps, in one degree of freedom), 
the expressions obtained are too complicated to be of real practical use. Even in one degree of 
freedom, the expressions up to fourth order are already very complicated. We therefore consider 
an iterative procedure for composing the generating functions, which simplifies the expressions 
required. 

2.4 Generating function composition: iteration method 

As an alternative to the induction method, we can attempt to solve the equations involved in gen- 
erating function composition using iteration. The formulae required may be simpler to imple- 
ment in practice. Suppose that we have two mixed variable generating functions, Fa{x\,p X 2) and 
Fb(x2,Px3), which are known. Fa gives the transformation from Oi to 02'. 

dF A 

x 2 = , (2.37) 



dx\ ' 



dF A 

Pxi = ^, (2.38) 



and F B gives the transformation from 02 to CT3 : 



x 3 = 3 , (2.39) 

a Pxi 

Px2 = ^- (2.40) 

0x2 

We wish to find a generating function Fc(x\,p x s) which gives the transformation from Oi to 03: 

x 3 = (2.41) 

aPxi 

dFr 

Pxi = -5-^. (2.42) 
dx\ 



We have seen in Section 2.3 that if the generating functions are expressed as polynomials in the 



appropriate variables, then the problem can be reduced to one of finding solutions to a set of poly- 



nomial equations of the form ( |2.2q ). We can attempt to find a solution iteratively, using Newton's 
method, as follows. First, we construct the matrix /: 

/ d3*i d» x \ 



ox 2 dp x2 
dx 2 dp x2 



(2.43) 

Then, if we make an initial estimate of the solution, for example: 

X (1) =r), (2.44) 



(where X is given by Eq. (2.28)) then we can construct improved estimates: 



X(„+i) = X (n) - 7 (n j J J j (2.45) 
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starting from n = I. For Eq. ( 2.26| ), the solution needs to be performed using an algebraic code, 
since the polynomials are expressed in terms of more variables than just those for which a solution 
is sought. However, we find in practice (e.g. for the illustration presented in Section |J), that the 
iteration can be carried out quickly and converges rapidly, even when working to high order in the 
variables. 

This approach to solving a system of nonlinear equations can also be applied to the Runge- 



Kutta integrator described in Section 2.2. If the Hamiltonian is expressed as a polynomial in the 
dynamical variables (note that this means expanding the square root in Eq. (p77|), but this can be 



done to higher order than is done for the paraxial approximation), then Eq. ( [2.1 6| ) represents 2n 
polynomial equations in the 2n "intermediate" variables {x' ,p' x ). These may be solved iteratively 
using Newton's method (in an algebraic code), to express the intermediate variables as power series 
in the initial variables {x\,p x i). Then, it is only necessary to make a substitution of p xl for p x \, 



using Eq. ( [2.1 5b| ), to produce an expression for the mixed-variable generating function ( J2. 19| ) for a 
single integration step. 



3. Example: quadrupole with octupole component and fringe field 

As an example, we illustrate the construction of the map, in the form of a mixed- variable generating 
function of the second kind, for a quadrupole with strong octupole component. The axis of the 
magnet (on which the field is zero) defines the reference trajectory. We allow both the quadrupole 
and octupole components to vary with position along the axis of the magnet, so as to represent 
some fringe field region at the entrance and exit of the magnet. An appropriate representation of 



the field, including the longitudinal variation, is given by Dragt's generalised gradients [13]. The 



appropriate expressions are given in reference [ 27 ] ; for convenient reference, they are reproduced 
in Appendix [B| In the notation used there, our field is defined by: 



C 2 {s) 



1 jfei 
2Bp 



sin (k s s) 



1 fc3 ■ 2/, \ 
4!^ Sm (M ' 



(3.1) 
(3.2) 



with k\ = — 10 m~ 2 , kj = 6 x 10 4 m~ 4 , and k s = 10m _1 . The total length of the magnet is defined 
to be L = %jk s . The expressions for the vector potential involve infinite summations: in practice, 
these must be truncated at some order in the co-ordinates. Although the fields corresonding to 
the truncated potentials strictly speaking do not satisfy Maxwell's equations, the error is of order 
higher than the truncation; and if the truncation is at significantly high order, the transfer map, itself 
computed only up to some order, is not affected. Note that the expressions for the transformations 
in the symplectic integrator involve the vector potential normalised by Po/q: therefore, it is not 
necessary to specify either the particle charge or the reference momentum. The normalised vertical 
field along the line x = — 1 mm, y = from s = to s = L is shown in figure |]. 



Using the generalised gradients ( p.l[ ) and (p.2D, we compute the mixed-variable generating 
function for the transfer map from s = to s = L, using the fourth-order Runge-Kutta integra- 



tor (Section 2.2), and the iterative method for composing the mixed- variable generating functions 



(Section 2.4). The Hamiltonian is given by Eq. (2.7), with the square root expanded to sixth order 
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0.05 0.1 0.15 0.2 0.25 0.3 
s(m) 



Figure 1. Vertical normalised field component in the example magnet, along the line x = — 1 mm, y = 
from s = to s = L. 



in the dynamical variables. Terms up to sixth order in the co-ordinates are retained in the expres- 
sion for the vector potential. We retain terms to fourteenth order in the dynamical variables when 
composing the generating function, to reduce feed-down errors. A step size of Aa = L/1024 was 
used in the integration. All computations were performed in Mathematica. 

The transfer functions (final co-ordinate and momentum, as functions of initial co-ordinate 
and momentum) are shown in figure g For comparison, we also computed the transfer functions 



using direct numerical integration of Hamilton's equations, using the Hamiltonian ( |2.7[ ) (without 
any approximation for the square root). The residuals between the transfer functions obtained from 
the mixed-variable generating function, and the transfer functions obtained from the numerical 
integration are also shown in figure ||. 

The residuals appear small over the range of variables plotted. To investigate the residuals 
more closely, we look specifically at the value of p x (L) as a function of x(0), for p x (0) = 0: we 
refer to this as the 1-D (one-dimensional) transfer function. A plot of the 1-D transfer function is 
shown in figure f| (note that there is some third-order curvature evident, from the octupole compo- 
nent in the quadrupole). Also shown in figure g is the difference (the residual) between the 1-D 
transfer functions obtained from the mixed-variable generating function and from direct numerical 
integration. The residual can be fitted (apart from some fluctuations, which may be due to limits 
on numerical precision) using a polynomical with only a single term, in x 1 . This indicates that the 
residual is dominated by the seventh-order terms in the map. 



We can also compare the coefficients in the 1-D transfer function (shown in figure g, left): 

5 

I 

m=\ 



Px (L)=Y J h m x(0) m . (3.3) 



The coefficients h m can be obtained using various techniques; a comparison of the results from 
some methods are shown in table [lj First, we obtain the h m directly from the generating function, 
by solving: 

dF 

P X l=0=j—, (3.4) 

ax\ 



10 " 10 



x(0) (mm) 

10 " 10 




Figure 2. Transfer functions in the example magnet. Top: horizontal co-ordinate (left) and normalised 
momentum (right) at s = L as functions of co-ordinate and momentum at s = 0. Bottom: residuals for the 
final horizontal co-ordinate (left) and normalised momentum (right), between the values calculated using the 
mixed-variable generating function, and a purely numerical integration through the field. 



where F is the generating function for a map through the entire magnet; the solution may be 
expanded in a power series in x\, which gives the h m directly. Second, we can obtain the h m from 
a polynomial fit to the 1-D transfer function obtained by numerical integration. Finally, estimated 
values for the h m may be obtained by integrating the appropriate normalised multipole component 
along the reference trajectory: 

h m k, — - / k m ds. (3.5) 
ml J 

The results from the first two methods are in good agreement up to the fourth-order Q14), and in rea- 
sonable agreement for The values obtained using the third method (integrating the normalised 



multipole component along the reference trajectory, Eq. (3.5)) are only approximate: effects asso- 



ciated with the fringe fields, and the non-zero length of the magnet, are not properly taken into 




x(0) (mm) x <°) (mm) 



Figure 3. Left: normalised horizontal momentum at s — L, as a function of horizontal co-ordinate at s = 0, 
for p x (0) = 0. Right: residual in final normalised horizontal momentum between values calculated from the 
mixed-variable generating function and from direct numerical integration. The black line shows the residual, 
the red line shows a fit to the residual, using a function of the form Ap x (L) x(Q) 7 . 



Table 1. Coefficients in the transfer function, derived from the mixed-variable generating function (MVGF), 
and by a polynomial fit to the tracking results using direct numerical integration (NI). Also shown for com- 
parison are the values that would be expected from the field integrals (FI). 



m 


h 


h 


h m 




MVGF 


NI 


FI 


1 


1.65228 m- 1 


1.65226 m- 1 


1.5708 m- 1 


2 


-1.00075 xl(T 13 rrr 2 


-1.00075xl0- 13 m- 2 





3 


-1933.15 m- 3 


-1930.82 m" 3 


-1570.8 m- 3 


4 


2.21872xl0- 9 m- 4 


2.21798 xl0- 9 m- 4 





5 


3.84174xl0 5 m- 5 


3.30479 xlO 5 m- 5 






account. 

Finally, we can see the effect the paraxial approximation would have if applied to this exam- 
ple, by comparing the residuals between the results of the direct numerical integration using the 



Hamiltonian (|2.7j), and the Hamiltonian in the paraxial approximation ( |2.10[ ). The difference be- 
tween the two is small in this case, but visible if we compare the difference in the final co-ordinate, 
as a function of initial co-ordinate and momentum: this is shown in figure ||. The residuals are 
an order of magnitude larger than those between the direct numerical integration for the "exact" 
Hamiltonian, and the integration (with mixed-variable generating functions) for the Hamiltonian 
expanded to sixth order in the dynamical variables. This suggests that using mixed-variable gener- 
ating functions, we can readily construct an exactly symplectic high-order map to good accuracy, 
without needing to make the paraxial approximation in the Hamiltonian for the system. 

A. Proof of the inductive formula for solving a system of polynomial equations 



We wish to show that (2.33) is a solution to (2.29). Assuming that Ti is invertible, we multiply 
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Figure 4. Residuals in the final co-ordinate (left) and final momentum (right) as a function of initial co- 
ordinate and momentum. The residuals are calculated as the difference between the results found from 
the mixed-variable generating function for the "exact" Hamiltonian, and the numerical integration of the 
Hamiltonian in the paraxial approximation. 



Eq. ( g29j ) by -f j" 1 to give: 



2V 

I 

r=2 



^T r (X,...,X)-X + T = O(e N+1 ). 



Let S„ be the left hand side of (A..1): 



(A.l) 



N 



s ;i = £ty(x,...,x)-x+t . 

r=2 



Substituting for X from Eq. ( |2.32| ) gives: 



(A.2) 



S re = £ T r ( £ c m , . . . , £ <VJ - ^ c,- + T . 

r=2 m\ = \ m T =\ r=\ 



Since ci = To, we find: 



N I N N \ N 

Sn = £T r £ c mi ,..., £ c Wr -£c r + 0(e' 

r=2 \mi=l m r =l / r=2 

N N N N 

= XI-LT r (c mi ,...,c m ,)-^ r + on 

r=2mi = l m r =\ r=2 

= IE I T r (c mi ,...,c mr )-f c, + 0( £ w+1 ) 

r=2k=2 m\,...^n^=\ k=2 
m,H Ym k =k 

~ l(" c * + I I T r (c mi ,...,c m j)+0(e iv + 1 ) 

&=2 \ r=2 mi,...,mj=l / 



(A.3) 



(A.4) 
(A.5) 
(A.6) 

(A.7) 



mi H hm,t=& 
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For r > k the set {(mi,..., m r ) m, > 1 and m\ + . . . +m r = k} = 0, hence: 

N / k 



S„ = £(-c fe +£ £ (c /l ,... J c/,)J+0( 

k=2 \ r=2 mi,...,mt=l / 



g^ 1 ). (A.8) 



mi H h»!j=/: 



Substituting for from Eq. (2.33) gives 



S n = 0(e N+1 ), (A.9) 



which is Eq. (2.29) 



B. Generalised gradients 

In Cartesian variables, a vector potential in the Coulomb gauge can be expressed as [27]: 



< co oo ■ 

1 CX3 OO I 

^ = i£3c+»r 1 ^(-i)' w , (< ;-, +1) , <ff H1 w^ + >')'. (b.2) 

oo oo i 



where: 



C« = (B.4) 

Maxwell's equations in free space are satisfied for any differentiable functions C m (z), which are 
called the generalised gradients. The conventional multipole expansion is obtained in the case that 
the C,„(z) are independent of z; then, Ci represents the dipole field component, C2 the quadrupole 
field component, and so on. 

The symplectic integrator described in Section^ uses a vector potential in the gauge: 

A y = 0. (B.5) 

This is readily obtained from the above expressions using the gauge transformation: 

A'=A-Vy, (B.6) 

where: 



Y = J A y dy. 



(B.7) 
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